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1  Introduction 


Microcalcification  detection  is  the  hallmark  of  mammography  as  a  breast  cancer  screening 
modality.  For  technical  reasons,  ultrasonic  detection  of  all  mammographically-visible  micro- 
calcihcations  has  been  problematic.  In  clinical  ultrasound,  high  frequencies  must  be  used 
to  resolve  microcalcihcations  below  200  micrometers.  Unfortunately,  ultrasonics  above  10 
MHz  suffer  from  appreciable  attenuation  in  soft  tissues,  and  depth  of  penetration  is  limited. 
Transmission  diffraction  tomography,  while  well-suited  for  the  geometry  of  the  breast,  is 
inherently  insensitive  to  scattering  caused  by  small,  hard  inhomogeneities.  A  more  general 
form  of  acoustic  inverse  scattering  is  therefore  needed  for  microcalcihcation  detection  and 
localization  by  ultrasound.  We  hnd  rationale  in  the  advanced  scalar  inverse  scattering  the¬ 
ory  developed  by  Colton,  Kirsch,  and  others  in  the  RADAR  community  that  can  determine 
the  shape  of  scatterers  with  size  on  the  order  of  the  wavelength.  In  addition  to  size  and 
number,  the  morphology  of  breast  microcalcihcations  is  an  important  diagnostic  indicator. 
Our  hypothesis  is  that  the  linear  sampling  method  (LS),  when  augmented  with  a  method 
for  estimating  the  inhomogeneous  Green’s  function  for  wave  propagation  in  the  breast,  can 
be  translated  to  an  acoustic  imaging  system  to  detect,  localize,  and  characterize  microcalci¬ 
hcations  in  breast  phantoms  using  data  from  the  scattering  measurements  in  a  tomographic 
geometry. 

2  Body 

The  goal  of  this  research  endeavor  is  to  develop  a  bistatic  ultrasound  imaging  method  that 
specihcally  targets  breast  microcalcihcations.  By  bistatic  imaging,  we  mean  that  receiver 
and  transmitter  can  be  separated  in  space.  Since  there  are  several  commercial  breast  acoustic 
tomography  systems  currently  undergoing  FDA  trials,  we  believe  that  it  is  the  appropriate 
time  to  apply  state-of-the-art  methods  from  optimum  array  processing  and  inverse  scattering 
to  this  important  biomedical  imaging  problem. 

In  year  2  of  this  research,  the  majority  of  our  ehorts  have  focused  on  development  of  a 
microcalcihcation  detection  algorithm  that  incorporated  a  priori  information  on  soft  tissue 
inhomogeneities  and  experimentation  in  phantoms  to  demonstrate  the  potential  of  the  pro¬ 
posed  method.  In  Task  6  of  the  original  Statement  of  Work,  we  have  continued  our  research 
into  developing  an  elliptical  Radon-based  approach  to  imaging  the  hetergeneous  background 
Green’s  function.  As  reported  in  year  1,  the  published  work  of  a  French  group  [1]  is  clearly 
incorrect  since  it  does  not  reproduce  the  known  results  for  spherical  Radon  transform  in  the 
limiting  case.  The  applied  mathematicians  supported  under  this  Synergistic  Idea  award  have 
made  signihcant  progress  in  analyzing  this  model,  and  a  paper  on  this  work  is  expected  to  be 
submitted  in  the  proposed  no-cost  extension  period.  In  summary,  we  have  found  some  suc¬ 
cess  in  diagonalizing  the  elliptical  Radon  transform  in  a  manner  that  produces  a  series-based 
inversion.  The  coefficients  in  this  reconstruction  algorithm  are  given  by  integral  equations 
that  can  be  solved  numerically.  Gurrent  efforts  by  collaborators  at  UT  Arlington  at  in  the 
direction  of  efficiently  implementing  this  inversion  as  a  numerical  method  that  can  be  trans- 
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lated  to  the  prototype  system  at  UT  Southwestern.  For  a  further  review  by  collaborators  on 
the  sub-award,  please  see  the  attached  Appendix  A. 

The  second  major  effort  for  year  2  of  this  research  project  was  imaging  experiments  with 
phantoms  and  simulated  microcalcihcations.  Six  high  signal-to-noise  multi-static  response 
matrix  data  sets  were  acquired.  These  data  supplement  the  scattering  data  from  wires  ob¬ 
tained  in  year  1  (so-called  point  scatterers),  and  inclnde  metal  scatterers  with  both  circular, 
hexagonal,  and  star  cross  sections.  The  strong  scatterers  were  embedded  in  a  7  cm  diameter 
gelatin  cylinder  with  known  aconstic  properties.  Image  reconstrnction  and  data  analysis  for 
these  data  sets  is  on-going  and  will  be  snbmitted  for  pnblication  in  the  proposed  no-cost  ex¬ 
tension  period.  In  addition,  since  Jnne  2009  we  have  been  investigating  more  sophisticated 
phantoms  for  validation.  In  one  case,  water-filled  voids  are  introdnced  into  the  gelatin  to 
prodnce  a  heterogeneons  backgronnd.  Flnctnations  in  the  data  dne  to  these  hetergeneities 
can  be  similar  in  appearance  to  the  signals  dne  to  the  strong  scatterer.  In  a  parallel  de¬ 
velopment  to  snpport  Task  7  (statistical  characterization  of  speckle  noise  and  snppression 
methods),  we  have  experimented  with  phantoms  with  clntter  to  prodnce  the  type  of  speckle 
observed  in  real  soft  tissne.  We  have  adopted  a  method  from  the  literatnre[2],  and  have  con- 
hrmed  the  B-mode  appearance  of  this  phantom  using  a  VisnalSonics  Vevo700  small  animal 
nltrasonnd  system.  In  the  proposed  no-cost  extension  period,  additional  data  acqnisition 
experiments  will  be  performed  as  a  part  of  Task  7  where  hard  2D  scatterers  will  be  inclnded 
in  phantom.  We  anticipate  completion  of  these  experiments  in  the  next  year  with  the  goal  of 
nnderstanding  the  effect  of  speckle  noise  on  the  the  robnstness  of  the  linear  sampling  method 
for  estimating  microcalcihcation  morphology. 

For  Task  13,  a  peer-reviewed  conference  proceeding  pnblication [3]  was  printed  in  early 
2009  (see  Appendix  B).  This  work  was  in  snpport  of  Tasks  6  and  12,  and  was  made  pos¬ 
sible  by  the  snpport  of  a  Biomedical  Engineering  graduate  student  through  this  research 
project.  In  addition,  a  collaborator  snpported  nnder  the  snb-award  pnblished  a  peer-reviewed 
book  chapter  on  aconstic  tomography  in  breast  imaging.  His  snpport  nnder  this  research 
project  was  critical  to  the  completion  of  this,  and  it  is  noted  in  the  acknowledgements  of 
the  pnblication [4].  At  present,  the  BME  gradnate  stndent  and  PI  at  UT  Sonthwestern  are 
preparing  an  expanded  pnblication  on  bistatic  nltrasonnd  imaging  of  scatterers  in  the  breast 
for  an  IEEE  transactions,  and  a  snpported  gradnate  stndent  in  Mathematics  and  collabo¬ 
rators  at  UT  Arlington  anticipate  at  least  2  snbmitted  pnblications  in  the  proposed  no-cost 
extension  period. 

In  the  next  year,  we  anticipate  completion  of  the  image  reconstruction  and  data  analysis 
tasks,  along  with  some  additional  experiments  to  test  the  robustness  of  the  described  methods 
to  snb-wavelength  scatter.  In  the  proposed  no-cost  extension  period,  snpport  will  continnte 
for  gradnate  stndents  at  UT  Sonthwestern  and  UT  Arlington,  while  the  PI  and  collaborators 
will  continne  the  Statement  of  Work  research  throngh  departmental  cost  sharing. 
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3  Key  Research  Accomplishments 

•  Data  collection  of  six  2D  phantoms  with  materials  mimicking  breast  tissues  and  mi- 
crocalhcations. 

•  Development  of  elliptical  Radon  transform  model  for  bistatic  estimation  of  background 
Green’s  function. 

•  Identihcation  of  fundamental  shortcoming  in  the  breast  ultrasound  tomography  litera¬ 
ture  as  it  relates  to  bistatic  imaging,  and  the  development  of  an  iterative  approach  for 
estimating  the  inhomogeneous  background  in  the  inverse  scattering  problem. 

4  Reportable  Outcomes 

•  A  publication  on  ultrasound  tomography  using  the  elliptical  Radon  transform  model[3]. 

•  A  review  article  on  tomographic  imaging  using  spherical  Radon  transform  models [4]. 

•  Two  collaborators  (M  Lewis  and  G  Ambartsoumian)  participated  in  the  2008  NSF/GBMS 
Rice  Workshop  on  Imaging  in  Random  Media,  where  image  reconstructions  methods 
of  the  types  of  interest  in  this  research  study  were  the  focus. 

5  Conclusion 

We  have  maintained  a  reasonable  focus  with  our  original  Statement  of  Work.  We  self-identify 
delayed  progress  on  Task  12  (Image  reconstruction  and  data  analysis)  which  will  be  remedied 
by  ongoing  work  in  the  no-cost  extension  period.  In  addition,  further  experiments  and  data 
acquisition  using  more  complicated  phantoms  will  faciliate  extensions  of  Tasks  7  and  11.  We 
anticipate  that  successful  completion  of  Task  12  will  produce  multiple  additional  publica¬ 
tions  that  will  support  translation  of  these  image  reconstruction  methods  to  a  commerical 
breast  ultrasound  tomography  system.  In  the  coming  year,  we  will  actively  solicit  research 
collaborations  in  this  area. 
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Appendix  A  -  Research  related  to  Task  6 
performed  under  UT  Arlington  subaward 


17  September  2009 


We  consider  a  simple  model  of  near-field  ultrasound  tomography,  where 
a  reflecting  medium  of  interest  is  excited  by  spherical  waves  emitted  from  a 
transducer,  and  the  backscattered  echoes  are  registered  by  a  receiver  [11,  12], 
Assuming  an  ideal  propagative  medium,  where  the  speed  of  sound  is  constant, 
the  signal  registered  at  any  given  moment  by  the  receiver  is  generated  by 
reflections  from  all  those  points  for  which  the  sum  of  their  distances  to  the 
emitter  and  the  receiver  is  constant  (depending  on  time  and  sound  speed). 
In  other  words,  those  points  are  located  on  confocal  ellipses  with  foci  at  the 
emitter  and  receiver  locations,  and  the  problem  of  image  reconstruction  in 
2D  boils  down  to  the  inversion  of  a  transform  integrating  functions  along 
such  ellipses.  In  3D  one  should  consider  a  transform  integrating  along  the 
family  of  surfaces  of  revolution  (obtained  by  rotating  the  2D  ellipse  around 
its  main  axis). 

Definition  1.  The  elliptical  Radon  transform  (ERT)  of  f{x),  x  E  'Mf  is 
defined  as 

Rf{Pe,Pr,r)=  j  f{x)dl{x), 

\x-pe\  +  \x-pr\=r 

where  dl{x)  is  the  arc  length  of  the  ellipse  \x  —  Pe\  +  \x  —  Pr\  =  r . 

It  is  easy  to  see  that  the  inversion  of  this  transform  is  an  overdetermined 
problem  (recovering  a  function  of  2  variables  from  a  function  of  5  variables). 
Hence,  one  should  expect  to  be  able  to  reconstruct  /  from  restrictions  of  Rf 
to  2— dimensional  families  of  ellipses.  From  physical  considerations  we  hrst 
restrict  the  locations  of  foci  Pe  and  pr  (the  acquisition  geometry)  to  a  curve  S 


in  (denoting  the  transform  by  Rsf{pe,Pr,r)  =  RfiPe:Pr,'>')\j,^^  Pr^s)-  We 
consider  the  lines  and  circles  for  S  in  this  project.  Notice,  that  Rsf{pe,Pr,  f) 
is  a  three-parameter  family,  while  f{x)  is  a  function  of  two  variables.  There 
are  several  ways  to  reduce  the  dimension  of  the  family  of  Rsf{pe,Pr:  f')-  Mo¬ 
tivated  by  tomographic  models  (e.g.  [16,  3,  11,  12])  we  will  assume,  that  the 
distance  between  the  emitter  and  the  receiver  stays  constant.  By  matching 
the  dimensions  this  way  we  further  denote  the  ERT  by  Rsf{p,  r),  where  S  is 
either  a  circle  or  a  line,  and  p  is  the  midpoint  of  the  interval  joining  the  foci 
of  the  ellipse. 

Notice,  that  in  the  degenerate  case  when  the  emitter  coincides  with  the 
receiver  the  ellipses  become  spheres  and  the  ERT  becomes  a  spherical  (or 
circular)  Radon  transform  (SRT).  The  latter  has  been  extensively  studied  in 
the  past  for  various  applications  of  imaging  and  applied  mathematics  (e.g. 
[1,  2,  4,  8,  9,  10,  13,  15]  and  the  references  therein).  In  case  of  SRT  various 
exact  inversion  formulae  have  been  discovered  for  a  limited  number  of  acqui¬ 
sition  geometries.  These  results  can  be  divided  into  two  categories:  closed 
backprojection  type  inversion  formulae,  and  expansions  into  series  (usually 
involving  some  special  functions).  One  of  the  goals  of  our  project  has  been 
obtaining  extensions  of  these  results  to  the  case  of  ERT,  and  the  discov¬ 
ery  of  some  new  inversion  algorithms.  While  the  work  in  this  direction  is 
continuing,  we  can  report  certain  progress  here. 

We  developed  an  inversion  algorithm  using  the  algebraic  reconstruction 
techniques  based  on  the  Kaczmarz  method,  which  can  work  in  any  acquisition 
geometry.  It  also  provides  an  effective  mechanism  for  incorporation  of  certain 
side  constrains  to  the  reconstruction  process,  which  is  an  extremely  important 
tool  for  stabilizing  the  inversion  in  limited  data  problems.  The  developed 
algorithms  have  been  numerically  implemented  in  Matlab,  and  are  currently 
being  tested  with  both  synthetic  and  experimental  data.  Two  publications 
describing  the  developed  methods  and  the  obtained  results  are  in  preparation. 
The  mathematical  description  of  the  inversion  instabilities  in  limited  data 
problems  has  been  covered  in  [17]. 

We  have  also  made  certain  progress  in  generalizing  the  inversion  formu¬ 
las  for  SRT  using  series  expansion  to  the  case  of  ERT.  Most  of  the  inver¬ 
sion  formulas  using  series  expansions  for  SRT  are  based  on  the  fact,  that 
this  transform  is  either  shift-invariant,  or  can  become  shift-invariant  after 
a  smooth  change  of  variables.  For  example,  the  result  of  [13]  is  derived  as 
follows.  Let  Rs  be  the  2D  spherical  Radon  transform  on  the  plane  that  inte¬ 
grates  functions  compactly  supported  inside  the  unit  disk  D  over  all  circles 
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\x  —  p\  =  p  with  centers  p  located  on  the  unit  circle  S.  Since  this  transform 
commutes  with  rotations  about  the  origin,  the  Fourier  series  expansion  with 
respect  to  the  polar  angle  partially  diagonalizes  the  operator,  and  thus  the 
n-th  Fourier  coefficient  gn{,p)  oi  g  =  Rsf  will  depend  on  the  n-th  coefficient 
fn  of  the  original  /  only.  Hence,  the  problem  of  reconstructing  /  from  g, 
can  be  reformulated  as  a  problem  of  reconstructing  /„  from  gn.  The  lat¬ 
ter  requires  solving  and  integral  equation,  which  in  [13]  is  done  using  the 
Hankel  transform.  In  [14]  the  SRT  is  inverted  in  the  case  of  a  linear  aper¬ 
ture  by  utilizing  the  shift  invariance  of  the  transform  in  the  direction  of  the 
data  acquisition  line,  and  obtaining  shift  invariance  in  the  orthogonal  direc¬ 
tion  by  a  smooth  change  of  variable.  Then  application  of  the  2-dimensional 
Fourier  transform  diagonalizes  the  operator,  which  enables  the  inversion.  In 
case  of  a  non- degenerate  ERT  in  linear  acquisition  geometry,  we  have  proved 
that  there  does  not  exist  any  smooth  change  of  coordinates,  which  would 
make  the  transform  shift-invariant  in  the  second  variable,  hence  one  should 
not  expect  inversion  procedures  using  Fourier  techniques  in  this  case.  In 
the  case  of  spherical  acquisition  geometry  the  elliptic  transform  R  also  com¬ 
mutes  with  rotations  about  the  origin,  hence  the  inversion  method  using  the 
Fourier  series  expansion  may  allow  generalization  to  ERT.  We  have  reduced 
the  problem  of  recovering  /  from  g  to  the  problem  of  recovering  /„  from  gn- 
The  latter  is  a  non-trivial  integral  equation,  and  we  currently  work  on  its 
solution. 
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ABSTRACT 

This  work  is  concerned  with  the  numerical  implementation  of  a  reconstmction  algorithm  developed  to  recover  a  function 
from  its  spherical  means  over  spheres  centered  on  a  circle.  The  algorithm  is  experimentally  verified  by  simulations  using 
numerical  phantoms.  In  the  scheme  of  tomography,  acoustic  waves  are  generated  by  illuminating  an  object  with  a  short 
burst  of  radio-frequency  waves.  In  applications,  like  breast  cancer  imaging,  which  use  modalities  like  photo-acoustic 
tomography  (PAT)  that  model  the  acoustic  pressures  as  spherical  means,  data  are  measured  on  the  detectors  located  in  a 
circle  surrounding  the  object.  This  is  then  used  to  reconstruct  the  absorption  density  inside  the  object.  In  contrast, 
applications  like  bore  hole  tomography  and  improved  Intravascular  Ultra  Sound  (IVUS)  imaging  for  prostate  cancer, 
which  use  modalities  like  Radial  Reflection  Diffraction  Tomography  (RRDT),  a  ring  of  detectors  placed  exterior  to  the 
object,  collect  the  acoustic  waves  as  back-scattered  field.  This  work  uses  a  single  algorithm  to  reconstruct  functions  from 
data  collected  using  these  two  different  techniques  -  one,  by  placing  the  object  inside  the  ring  of  detectors,  and  the  other, 
by  placing  the  object  exterior  to  the  ring  of  detectors.  The  algorithm  then  draws  a  comparison  between  the  two 
reconstructions.  The  case  of  bistatic  ultrasound  imaging,  where  the  elliptical  Radon  transform  is  appropriate,  is  also 
discussed. 

Keywords:  Ultrasound  Tomography,  Spherical  Radon,  Elliptical  Radon,  Diffraction  Tomography 


1.  INTRODUCTION 

In  imaging  techniques  like  photo-acoustic  tomography  biological  tissues  [1]  are  irradiated  with  short  pulse  waves  of 
energy,  and  ultrasound  waves  are  subsequently  generated.  In  medical  modalities  such  as  breast  imaging  [2]  these  waves 
are  measured  by  sensors  on  the  exterior  surface  of  the  tissues.  Under  the  approximation  that  the  speed  of  sound,  c,  is 
constant  throughout  the  body,  in  a  monostatic  setup  (where  a  single  transducer  acts  as  both  the  source  and  detector),  the 
measured  ultrasound  waves  are  modeled  as  circular  means  of  their  reflectivity  obtained  over  concentric  circles  centered 
at  the  measured  point  (fig  1).  Formulas  of  the  filtered  back-projection  to  reconstruct  the  object  were  established  [3]  for  a 
function  supported  in  the  ball  in  even  dimensions.  In  this  paper,  we  consider  implementing  this  algorithm  to  reconstruct 
phantoms  with  no  noise. 

The  algorithm  is  also  implemented  to  synthetic  data  obtained  by  placing  the  ring  of  detectors  exterior  to  the  object.  In 
applications  like  Radial  Reflection  Diffraction  Tomography  (RRDT)  [4]  that  can  be  used  in  medical  modalities  such  as 
prostate  imaging  ,  the  object  is  external  to  transducer  geometry  (fig  2)  and  the  transducers  launch  a  primary  field  and 
collect  the  ultrasound  waves  as  backscattered  field.  In  such  multimonostatic  assemblies  (where  a  single  transducer  acts 
both  as  a  source  and  detector  at  multiple  locations)  the  reconstruction  techniques  involve  different  kind  of  wave 
modeling.  We  implement  the  formula  obtained  in  [3]  to  reconstruct  object  for  this  type  of  modality.  With  a  different 
method,  the  result  appears  to  be  similar  to  reconstructions  in  the  diffraction  tomography  technique  [4]. 

The  effect  of  the  algorithm  in  reconstructing  partial  data  is  also  analyzed.  Reconstruction  from  incomplete  data  when  the 
detectors  are  placed  along  a  semi-circle  surrounding  the  object  leads  to  interesting  observations  which  is  analytically 
accounted  for  in  [5]. 

Finally,  we  naively  implement  the  same  algorithm  by  slightly  changing  the  filter  to  reconstruct  data  obtained  in  a  bistatic 
setup  (one  transducer  acts  as  a  source  and  another  acts  as  the  detector).  This  method  uses  elliptical  Radon  transform  to 
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model  the  data.  There  is  no  filtered  back-projection  kind  of  algorithm  found  in  literature  that  reconstructs  elliptical 
Radon  transformed  data. 


Fig.  1  Object  placed  inside 
the  ring  of  detectors 


2.  GEOMETRICAL  SETUP 

In  near-field  ultrasound  tomography,  where  the  waves  are  spherical  and  the  sensors  are  point  sources/receivers,  various 
tomographic  geometries  are  possible.  Flere,  we  discuss  three  possible  2-D  geometrical  setups  in  the  case  where  the 
transducers  are  placed  along  a  circle.  An  equivalent  would  be  a  transmitter/receiver  being  rotated  about  a  fixed  center 
like  shown  in  fig  3. a.  At  each  point  on  the  circle  the  transducer  emits  a  wave  and  receives  the  scattered  wave  back,  an 
operation  called  as  “pitch/catch”  in  the  monostatic  case  [4].  In  the  2-D  case,  the  object  to  be  reconstructed  is  in  the  plane 
perpendicular  to  the  axis  of  rotation  of  the  transducer.  In  Radial  Reflection  Diffraction  Tomography  (RRDT)  algorithm, 
the  data  are  collected  as  backscattered  waves  at  various  angular  locations  and  the  inversion  is  performed  based  upon  a 
linear  scattering  model  [6]. 

In  this  paper  we  implement  an  algorithm  which  inverts  data  obtained  at  various  angular  locations  on  a  circle  as 
projections  which  are  spherical  means  (circle  in  2-D  with  the  transducer  location  as  the  center)  of  the  function  (object)  to 
be  reconstructed. 

2.1  Monostatic  setup 

As  discussed  above,  in  the  monostatic  setup  (fig.  4.a),  the  same  transducer  acts  as  both  the  transmitter  and  receiver. 
Assuming  the  speed  of  sound,  c,  constant  inside  the  object,  we  can  say  that  the  time,  t,  for  the  pulse  to  reach  the  medium 
and  time  the  echo  takes  to  reach  the  receiver  are  equal.  The  loci  of  all  points  in  the  medium  which  satisfy  the 
aforementioned  condition  will  lie  on  a  circle  with  ‘radius’ t.  Therefore  if  the  resultant  echoes  are  recorded  as  function  of 
time,  we  will  obtain  a  family  of  line  integrals  of  the  function  along  concentric  circles  centered  at  the  source  at  a 
particular  angle.  We  obtain  such  circular  integrals  at  multiple  angular  locations  along  a  circle.  Our  aim  is  to  reconstruct 
the  function  from  these  projections. 

Form  the  setup  explained  above,  we  can  parameterize  a  spherical  Radon  transform  with  two  parameters,  angle  0  and 
radius  r.  A  symbolic  statement  of  the  forward  equation  would  be  the  following: 


Fig. 2.  Object  placed  outside 
the  ring  of  detectors 


G(p,‘l’ )  =  l|R-r|/(r,0)  ds 


where  G(p,(|) )  is  the  line  integral  of  /  (r,0)  along  the  circle  of  radius  p,  whose  center  (the  source  transducer)  lies  on  the 
circumference  of  the  enclosing  circle  at  angle  (|) . 
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Fig. 3.  The  circular  integral  of  the 
function  is  obtained  over  a  circle 
centered  at  the  source  transducer,  r 
is  the  distance  from  the  origin  to 
the  scatterer  and  R  is  the  radius  of 
the  circle  along  which  the 
transducer  rotates. 


Fig. 4.  Various  geometrical  setups  for  a 
2-D  ultrasound  tomography  system 
where  the  transducers  are  placed  along 
a  circle  -a.  Monostatic  -  a  single 
transducer  acts  both  as  transmitter  and 
receiver,  b.  Bistatic  -  One  transducer 
acts  as  a  transmitter  and  other  acts  as 
receiver,  c.  Multistatic  -  One 
transducer  transmits  ultrasonic  waves 
and  many  transducers  act  as  receivers 
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2.2  Bi$tatic  setup 


As  shown  in  fig.  4.b,  in  a  bistatic  setup,  a  pair  of  transducers  moves  along  the  circle.  One  transducer  acts  as  a  transmitter 
and  the  other  acts  as  the  receiver.  If  we  consider  the  time  from  the  emitter  to  the  scatterer  is  tl  and  from  the  scatterer  to 
the  receiver,  t2,  and  assume  tl+t2  is  a  constant,  then  the  loci  of  all  scatterers  which  satisfy  this  condition  will  he  on  an 
ellipse.  In  this  case  the  recorded  resultant  echoes  at  the  receiver  as  a  function  of  time  will  be  line  integrals  of  the  function 
to  be  recovered  over  a  family  of  ellipses. 


Fig. 5.  The  line  integral  of  / 
in  the  case  of  bistatic  setup 
is  obtained  over  an 
elliptical  path.  An  ellipse  is 
defined  by  five  parameters. 


There  is  no  filtered  backprojection  (FBP)  kind  of  algorithm  in  literature  to  invert  an  elliptical  Radon  transformed  data.  A 
Fourier  transform  based  reconstruction  is  attempted  in  [7].  Flowever,  we  implement  the  same  filter  suggested  in  [3]  with 
a  slight  modification  to  invert  elliptical  Radon  transformed  data  in  a  FBP  fashion  for  a  bistatic  model. 

To  remark  on  the  parameterization  of  the  ellipse,  there  are  five  degrees  of  freedom.  But  to  make  the  model  practical  and 
physical  we  can  have  constraints  on  two  of  the  parameters  and  reduce  the  degrees  of  freedom  to  three.  The  problem  still 
remains  overdetermined. 

2.3  Multistatic  setup 


Multistatic  setup,  as  shown  in  fig.4.c,  has  one  source  and  more  than  one  receiver  and  this  setup  is  rotated  along  a  circle. 
Again  a  wave  based  tomographic  image  reconstruction  is  attempted  in  [4]. 


3.  FILTERED  BACKPROJECTION  ALGORITHM 


The  problem  of  reconstructing  the  function  from  the  circular  integrals  is  similar  to  the  classical  problem  of 
reconstructing  the  function  from  the  projections  in  x-ray  tomography.  In  fact,  the  circular  Radon  reconstruction  reduces 
to  classical  Radon  solution  in  the  limit  as  the  radius  of  the  as  the  radius  of  the  enclosing  circle  approaches  infinity. 
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We  were  primarily  inspired  by  the  FBP  algorithm  suggested  in  [3]  for  inverting  spherieal  Radon  transformed  data.  The 
paper  suggests  an  FBP  algorithm  with  linear  interpolation  in  two  dimensions.  The  logarithmie  filter  suggested  was 
dependent  on  the  radial  samples  of  the  aequired  data. 


The  algorithm  was  implemented  in  Matlab  and  for  all  (N+1)^  reeonstruetion  points,  N^+1  eenters  on  S  were  summed  up. 
If  there  are  pixels  and  if  N  is  approximately  equal  to  both  the  radial  and  angular  samples,  the  algorithm  requires  o(N^) 
operations. 


4.  METHODS  AND  RESULTS 

4.1  Case  1:  The  object  being  placed  inside  tbe  ring  of  detectors 

The  filtered-baekprojeetion  algorithm  suggested  in  [3]  was  applied  to  2-D  synthetie  data  obtained  from  performing 
forward  spherieal  Radon  transform  on  a  100  by  100  Shepp-Logan  phantom.  In  the  first  ease,  where  the  forward  Radon 
was  performed  with  the  objeet  inside  the  ring  of  deteetors  with  a  number  of  radial  samples  Nr=100  and  number  of 
angular  samples  Nphi=100.  The  sinogram  obtained  is  depieted  in  figure  6b  and  the  reeonstmeted  image  is  shown  in  6a. 


Fig. 6a.  Reeonstmeted  image  when  Fig.6b.  Corresponding  sinogram 

the  objeet  is  plaeed  inside 
ring  of  deteetors 


4.2  Case  2:  The  object  being  placed  outside  the  ring  of  detectors 


The  same  algorithm  was  applied  to  the  objeet  when  it  was  plaeed  outside  the  ring  of  deteetors  with  the  same  Nr  and 
Nphi.  The  sinogram  obtained  is  depieted  in  figure  7b  and  the  reeonstmeted  image  is  shown  in  7a. 
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Fig. 7a.  Reconstructed  image  when  Fig. 7b.  Corresponding  sinogram 

the  object  is  placed  outside 
ring  of  detectors 


4.3  Case  3:  Effect  of  incomplete  data 


We  tried  reconstmcting  the  object  by  placing  the  transducers  along  a  semi-circular  arc  at  the  bottom  of  the  object.  We 
used  a  different  phantom  to  understand  the  effect  of  algorithm  on  incomplete  data.  The  phantom  is  depicted  in  fig  8.a 
with  the  detector  locations  shown. 


Fig. 8. a.  Phantom  used  for 
reconstruction  from  partial 
data.  The  detectors  are  placed 
along  a  semicircular  arc  at  the 
bottom  of  the  phantom. 
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The  following  images  depict  the  sinogram  and  the  reconstructed  object: 
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Fig.  8.b  Sinogram  obtained  for  incomplete  data 


Fig.S.b  Reconstructed  image 


4.4  Case  4:  Applying  a  modified  filter  to  reconstruct  elliptical  Radon  transformed  data 


We  finally  implement  the  same  algorithm  with  a  slightly  modified  filter  to  reconstruct  elliptical  Radon  transformed  data. 
As  explained  in  section  2.2,  an  ellipse  has  five  degrees  of  freedom.  We  keep  the  distance  between  the  foci  constant 
(transmitter  and  receiver)  and  make  it  twice  the  dimension  of  the  Shepp-Logan  phantom  we  used  previously.  We  also 
make  the  center  of  the  foci  he  on  a  circle.  We  also  have  an  additional  constraint  on  the  rotation  of  the  axis  of  the  line 
connecting  the  foci.  We  keep  this  axis  perpendicular  to  the  radius  of  the  enclosing  circle. 

The  filter  was  modified  by  replacing  the  radial  factor  by  minor  axis  of  the  ellipse.  Projections  contained  samples  of 
minor  axes  at  every  angular  location.  We  used  the  Shepp-Logan  phantom  used  previously.  Fig  9. a  shows  the  sinogram 
obtained  and  9.b  shows  the  reconstructed  image. 


Fig. 9. a.  Sinogram  of  the  elliptical  Radon  Fig.9.b  The  reconstructed  image  with 

transformed  data  the  modified  filter 
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5.  CONCLUSIONS 


To  conclude,  we  set  out  to  implement  the  same  reconstruction  algorithm  for  the  same  mode  of  data  acquisition  with 
different  modeling  techniques.  In  the  scheme  of  photo-acoustic  tomography,  as  is  the  first  case  described  above,  the 
reconstruction  works  very  well,  whereas,  in  the  model  where  the  object  is  placed  outside  the  ring  of  detectors,  the 
algorithm  seems  to  reconstruct  only  the  artifacts. 

With  incomplete  data,  where  the  transducers  are  placed  in  a  semicircular  arc  around  the  object,  the  reconstruction  show 
similar  results  as  explained  in  [5].  We  observe  the  boundaries  disappear  in  the  case  of  incomplete  data.  Reconstruction 
from  elliptical  Radon  transform  data  gives  a  low  intensity  image  with  a  modified  fdter. 

There  is  no  existing  FBP  kind  of  algorithm  for  elliptical  Radon  transformed  data.  That  can  be  considered  for  a  future 
research  work.  The  problem  of  exterior  reconstruction  can  also  be  addressed  as  they  have  potential  impact  in  medical 
imaging  of  breast  and  prostate  cancer. 
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Abstract 

Photo-  and  thermo-acoustic  tomography  are  non¬ 
ionizing  imaging  techniques  that  reconstruct  internal 
photo  (thermo)  -acoustic  source  distributions  from  data 
acquired  as  projections  by  ultrasound  detectors  placed 
over  a  surface  tliat  encloses  the  object  under  study. 
They  combine  the  high  contrast  in  electromagnetic 
(EM)  absorption  between  healthy  and  cancerous  tissue 
with  the  high  resolution  of  ultrasound.  This  mode  of 
tomographic  imaging  has  recently  provoked  an  interest 
in  researchers  working  on  tlieoretical  aspects  of  image 
reconstruction.  Diffraction  tomography  methods  for  the 
same  geometry  exist  and  use  the  Bom  approximation 
for  the  reconstruction  techniques. 

We  investigate  a  Spherical  Radon  transform-based 
algorithm  developed  to  recover  a  function  from  its 
spherical  means  over  spheres  centered  on  a  circle.  The 
algorithm  is  numerically  implemented  and 
experimentally  verified  using  numerical  phantoms.  We 
implement  the  same  algorithm  to  reconstruct  data  from 
projections  obtained  by  placing  the  object  inside  the 
ring  of  detectors  and  placing  it  exterior  to  it.  Both  these 
arrangements  are  significant  as  they  represent  imaging 
techniques  like  Photo-Acoustic  tomography  and  Radial 
Reflection  Diffraction  Tomography  for  endo-rectal 
prostate  imaging. 


Rationale 

Anticipating  the  introduction  of  clinical  breast 
acoustic  tomography  systems  where  attempts  are  being 
made  to  solve  the  complete  non-linear  inverse  scattering 
problem,  we  are  considering  a  mathematical  model  that 
expresses  data  acquired  as  spherical  Radon  transformed 
projections.  The  idea  behind  this  monostatic  modality 
(the  same  sensor  is  used  as  both  transmitter  and 
receiver)  is  to  illuminate  an  object  by  a  short  burst  of 
radiofrequency  waves  which  cause  thennal  expansion 
and  subsequently  generates  an  acoustic  wave.  The 
existing  algorithms  in  both  photo-acoustic  and 
diffraction  tomography  successfully  reconstruct  the 
object  if  the  acoustic  waves  are  measured  on  the 
periphery  or  in  the  exterior  of  the  object.  Does  the  same 
algorithm  work  if  the  ring  of  detectors  is  placed  exterior 
to  the  object?  This  is  of  interest  because  some 
anatomical  constraints  like  in  the  case  of  prostate 
imaging  require  us  to  use  an  "exterior”  geometry. 
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Image  Reconstructions  for  Various  Geometrical  Configurations 


The  upper  half  of  the  quadrant 
shows  the  reconstruction  of  a 
101X101  grayscale  ‘Shepp- 
Logan  Phantom  when  the 
phantom  is  placed  inside  the 
ring  of  detectors. 


The  bottom  half  shows  the 
reconstruction  of 
The  same  phantom  when  the 
phantom  is  placed  outside  the 
ring  of  detectors 


Background 


Acoustic  imaging 
geometry 

Complex  amplitude  data 
for  each  receiver  Rj  and 
transmitter  Tj  is  organized 
in  multistatic  response 
matrix  Ky. 


f 


Forward  Problem 


■  The  Spherical  mean  transform  is  defined  by 
the  above  equation  (Finch,  et  ai.,  SIAM  J. 

Appi.  Math.  68(2007)  no.2,  392-412).  S^'^ 
denotes  the  area  of  the  unit  sphere  and  ds(0) 
denote  the  area  measure  on  the  sphere.  The 
aigorithm  was  performed  for  a  2-D  case  where 
n=2. 


The  algorithm  proposed  by 
Finch  et  at,  shows  good 
reconstruction  results  when 
the  object  is  placed  inside  the 
ring  of  detectors 


The  algorithm  seems 
to  reconstruct  only 
the  artefacts  when 
the  object  is  placed 
outside  the  ring  of 
detectors 


EtTect  of  Incomplete  data 


We  tried  reconstructing  the  object  by  placing  the 
transducers  along  a  semi-circular  arc  at  the  bottom  of 
the  object.  We  used  a  different  phantom  to  understand 
the  effect  of  algorithm  on  incomplete  data.  The 
phantom  and  the  reconstructed  image  are  depicted  in 
figures  below  with  the  detector  locations  shown. 


Potential  Impact 


•  Far  field  -  plane  wave  insonification,  asymptotic 
scattering 


Reconstruction  of  Elliptical  Radon 
Transformed  data  using  the  same  algorithm 


filtered  hack-projection  (FBP)  algorithm  (Finch, 
Haltmeier,  Rakesh,  SIAM  J.  Appl.  Math.  68(2007) 
no.2,  392-412)  with  linear  interpolation  in  two- 
dimensions  to  recover  a  function  from  its  spherical 
means.  Explicit  inversion  formulae  for  the 
spherical  radon  transform  is  proposed  by 
Kunyansky,  (Inverse  Problems.  23  (2007),  no.l, 
373-383). 

A  discrete  FBP  for  bistatic  modeling  involving 
elliptical  radon  transform  which  is  of  interest  to 
tomographic  imaging  is  not  yet  found  in  the 
literature. 


In  a  bistatic  setup  the  acquired  data  are  modeled  as  elliptical 
Radon  transform  projections.  If  we  assume  the  total  time  taken 
by  the  pulse  to  reach  the  scatterer  and  the  echo  to  reach  the 
receiver  a  constant,  then  the  loci  of  all  the  scatterers  which 
satisfy  this  condition  will  be  an  ellipse. 

We  changed  the  ‘radial’  factor  in  the  filter  in  the  original  FBP 
algorithm  to  use  the  minor  axis  of  the  ellipse  rather  than  the 
radius  of  the  projection. 


Ultrasound  is  not  currently  used  for  screening  of 
breasts  for  non-palpable  lesions;  but  with  proven 
sensitivity  and  specificity,  a  non-ionizing 
ultrasonic  method  has  the  potential  to  allow 
earlier  and  more  frequent  breast  screenings  in  at- 
risk  patients.  If  breast  compression  can  be 
avoided  by  the  imaging  technique,  then 
ultrasound  also  represents  a  less  stressful 
procedure  for  the  patient.  Breast  ultrasound  may 
be  the  only  available  technique  for  screening  in 
the  radiologically  dense  breasts  of  women  under 
40.  Advances  in  early  detection  have  already 
increased  survivability,  and  there  is  reason  to 
believe  that  improvements  in  breast  screening 
with  microcalcification-sensitive  ultrasonic 
imaging  will  continue  this  trend. 

There  are  no  existing  mathematical  model  in 
ultrasound  tomography  to  image  exterior  objects 
which  will  be  helpful  in  imaging  diseases  like 
prostate  cancer. 
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6.1  INTRODUCTION 

Thermoacoustic  tomography  (TAT)  employs  the  well-known 
[1-4]  correlation  between  electromagnetic  (EM)  absorption  of 
biological  tissue  and  its  physiological  and  pathological  proper¬ 
ties.  To  employ  this  contrast  mechanism,  a  biological  object 
is  irradiated  by  a  brief  EM  pulse,  and  the  resulting  thermoa¬ 
coustic  signals  from  the  tissue  are  collected  by  ultrasound 
transducers  to  map  the  distribution  of  the  radiation  absorption 
within  the  sample  (e.g..  Refs.  [5-9]).  TAT  thus  combines  the 
good  spatial  resolution  of  ultrasound  imaging  with  the  good 
contrast  in  EM  absorption. 


The  problem  we  address  in  this  chapter  is  of  limited  data 
(limited  view),  when  transducers  cannot  be  placed  around 
the  complete  object.  This  situation  is  very  common  in  TAT, 
for  instance  in  its  applications  to  breast  imaging,  where  only 
a  half  sphere,  rather  than  a  sphere  is  available  for  placement. 
The  question  that  arises  is  the  theoretical  and  practical  pos¬ 
sibility  of  reconstruction  from  limited  view  data.  Although, 
as  we  will  see,  mathematically  rigorous  uniqueness  results 
exist  that  guarantee  the  theoretical  possibility  of  recon¬ 
struction  from  limited  view  data,  the  practical  situation  is 
somewhat  different.  Namely,  some  features  of  the  object 
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to  be  imaged  (which  we  dub  “invisible”)  are  practically 
impossible  to  reconstruct,  because  they  suffer  from  a  man¬ 
datory  blurring.  Any  attempts  to  overcome  this  difficulty 
are  futile,  unless  some  prior  information  about  the  object  or 
missing  data  is  incorporated,  or  the  very  physical  set  up  of 
the  measurement  is  altered.  On  the  contrary,  the  “visible” 
details  are  easy  to  reconstruct  stably.  It  is  thus  important  to 
understand  these  limitations  due  to  limited  data.  Limited 
view  problems  have  been  studied  extensively  in  other,  more 
traditional,  types  of  tomography,  such  as  x-ray,  SPECT, 
diffraction  tomography,  and  reflectivity  tomography  (e.g.. 
Refs.  [10,11,13-18]).  The  TAT  situation  is  similar  and  has 
been  discussed  in  various  ways  in  Refs.  [19-27].  Although 
the  underlying  mathematical  technique  is  rather  involved, 
the  final  results  are  easy  to  understand  and  use.  The  goal 
of  this  chapter  is  to  review  these  results  and  to  provide  the 
relevant  references.  We  will  not  attempt  to  go  through  the 
rigorous  mathematical  technicalities,  but  rather  explain  the 
basic  ideas.  Correspondingly,  the  results  will  sometimes  be 
stated,  when  the  authors  feel  no  danger  of  misuse,  without 
some  necessary  technical  conditions  under  which  they  are 
proven.  References,  however,  will  be  given  where  the  rigor¬ 
ous  details  can  be  found. 

A  significant  limitation  of  what  is  described  in  this  chap¬ 
ter  is  that  the  speed  of  ultrasound  in  the  imaged  tissue  is 
assumed  to  be  constant.  The  case  of  a  variable  speed  could 
also  be  treated,  but  this  would  require  a  somewhat  different 
and  more  complex  consideration,  which  apparently  has  never 
been  completely  implemented. 

The  paper  is  structured  as  follows:  In  the  next  section, 
we  state  the  model  and  briefly  review  the  existing  unique¬ 
ness  results,  reconstruction  formulas,  and  procedures  for  the 
full  data  view  problem.  The  section  that  follows  contains  the 
description  of  how  the  limits  of  view  influence  the  unique¬ 
ness  and  reconstructions.  It  is  explained  how  one  can  deter¬ 
mine  which  interfaces  in  the  object  will  be  blurred  in  the 
reconstructed  image  due  to  the  limited  view.  We  also  intro¬ 
duce  in  this  section  necessary  simple  mathematical  notions. 
Then  the  next  section  presents  numerical  examples  (both 
from  synthetic  and  experimental  data)  that  illustrate  the  con¬ 
cepts.  All  mathematical  conclusions  are  equally  applicable  to 
photoacoustic  tomography. 

6.2  MATHEMATICAL  MODEL  AND 
RECONSTRUCTION  IN  TAT 

The  accepted  mathematical  model  of  TAT  involves  quite  a 
few  physical  constants.  As  it  happens,  their  presence  is  irrel¬ 
evant  for  understanding  the  concepts  and  using  the  limited 
view  results.  We  thus  present  here  a  simplified  model,  where 
all  constants  are  assumed  to  be  equal  to  one.  The  full-blown 
models  with  all  the  complications  (which  do  not  influence  the 
issue  we  discuss  here),  as  well  as  more  details  and  references, 
can  be  found  in  surveys  [28-30]  and  in  chapters  [31,32]  in 
this  volume. 

With  this  simplification,  the  model  of  TAT  is  described  by 
the  following  wave  equation  problem: 


Pa  =  L  ^  .r  G 

p(x,0)  =  f{x), 

p,(x,0)  =  0, 

P(y,t)=  g(y,t)  for  yeS,  t>0. 

Here  p{x,t)  is  the  pressure  at  the  point  x  and  time  t,  vXx)  is 
the  speed  of  the  ultrasound  propagation  in  the  tissue,  S  is  the 
observation  surface  where  the  transducers  are  placed,  g{y,t) 
is  the  measured  data,  i.e.,  the  value  of  the  pressure  at  time 
t  measured  at  the  transducer’s  location  y  eS,  and/(x)  is  the 
object  to  be  reconstructed.  As  we  have  already  mentioned, 
the  results  described  here  apply  to  the  case  of  a  constant 
sound  speed  only.  We  can  also  assume  that  v^=  1,  thus  arriv¬ 
ing  at  the  simpler  equations 

Pa  =  ^.P’ 

p(x,0)  =  f(x),  ^g2) 

p,(x,0)  =  0, 

_p(y,t)  =  g(y,t)  for  yeS,  t>0. 

In  this  (constant  speed)  case,  the  well-known  Kirchhoff- 
Poisson  formulas  (see,  e.g..  Refs.  [33,34])  allow  one  to  rep¬ 
resent  the  solution  pix,f)  in  terms  of  the  spherical  means  of 
the  function /(x) 


iRf)(y,t): 


Ajr 

ItoN 


f{y  +  t03)d(0 


for  yeS. 


(6.3) 


Thus,  knowledge  of  the  data  g(y,t)  is  equivalent  to  knowing 
the  integrals  of  the  unknown  function /over  all  spheres  cen¬ 
tered  at  transducers’  locations  yG  5.  One  immediately  notices 
a  similarity  with  the  standard  Radon  transforms  for  x-ray 
tomography  or  MRI,  where  the  integration  is  done  over  lines 
or  planes  rather  than  spheres.  And  indeed,  most  of  the  stan¬ 
dard  tomographic  results  and  techniques  find  their  analogs  in 
TAT,  albeit  many  aspects  become  much  more  involved  (see, 
e.g..  Refs.  [28-32,35-43]).  We  now  address  the  uniqueness  of 
reconstruction  and  reconstruction  procedures  in  TAT. 

The  reader  can  notice  that  in  all  feasible  applications,  func¬ 
tion /(x)  to  be  reconstructed  does  not  have  “infinite  tails”,  i.e., 
vanishes  outside  of  a  bounded  domain.  Moreover,  in  most 
practical  applications,  it  is  true  that  its  support  is  completely 
surrounded  by  the  observation  surface  S.  In  what  follows, 
this  will  be  assumed  (some  results  and  inversion  formulas 
might  fail  unless  this  is  satisfied  [30,31]). 

We  first  assume  the  knowledge  of  the  full  data  g{y,t)  for 
an  observation  surface  S,  which  is  a  sphere  surrounding  the 
support  of  the  image /(x).  The  case  of  limited-angle  data  will 
be  considered  in  the  next  section. 

A  similar  problem  to  the  one  we  have  just  described  in  3D 
can  also  be  considered  in  other  dimensions.  For  TAT,  only 
dimensions  2  and  3  are  relevant  (the  need  to  use  a  2D  prob¬ 
lem  arises,  for  instance,  when  one  uses  linear,  rather  than 
point-like,  detectors  [24,25]). 
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6.2.1  Uniqueness  oe  Reconstruction 

The  first  question  to  answer  is  whether  the  data  determines 
function  /  uniquely,  i.e.,  which  observation  surfaces  S  pro¬ 
vide  uniqueness  of  reconstruction.  In  2D,  this  was  resolved  in 
Ref  [35].  There  is  still  no  complete  solution  in  3D  (see  Refs. 
[15,28-32,35,38-41]  for  surveys  and  references).  However, 
from  the  practical  point  of  view,  there  is  no  problem.  Indeed, 
it  has  been  known  at  least  since  Ref  [44]  (see  also  the  ref¬ 
erences  just  mentioned  and  Ref.  [45])  that  if  the  surface  S 
is  closed,  e.g.,  a  sphere,  and  if  /  has  a  bounded  support  (not 
necessarily  surrounded  by  S),  then  there  is  uniqueness.* 
When  5  is  a  cylinder,  uniqueness  also  holds,  and  if  5  is  a 
plane,  uniqueness  holds  if/vanishes  on  one  side  of  the  plane 
(otherwise  an  odd  function  /  with  respect  to  the  plane  pro¬ 
vides  a  counterexample).  This  essentially  covers  any  feasible 
full-view  TAT  situation. 

Now,  when  uniqueness  is  established,  one  wonders  how  to 
actually  reconstruct  the  image/ 

6.2.2  Reconstruction  Formulas  and  Algorithms 

Although  the  area  had  experienced  a  slow  start,  there  has 
been  a  large  variety  of  inversion  formulas  and  algorithms 
developed  lately  (at  least  for  the  case  of  a  constant  sound 
speed  that  we  consider  here).  One  can  find  a  thorough  discus¬ 
sion  of  inversion  formulas  and  reconstruction  algorithms  in 
the  surveys  [31,32]  in  this  volume,  as  well  as  in  Refs.  [29-32]. 
We  will  just  say  a  few  general  words  about  these,  since  the 
details  of  a  particular  algorithm  do  not  seem  to  affect  the 
general  conclusion  that  we  will  make  later  about  the  limited 
view  problems  in  TAT. 

Explicit  inversions  based  on  Fourier  transform  techniques 
have  been  developed  in  the  case  of  a  planar  observation  sur¬ 
face  S  (see,  e.g..  Refs.  [16,46-50]). 

The  most  interesting  case  of  a  spherical  surface  S  was 
first  rigorously  treated  in  Refs.  [51,52],  where  the  rotational 
invariance  of  the  problem  was  used  to  expand  the  data  and 
the  image  into  Fourier  series  with  respect  to  the  angle.  This 
resulted  in  explicitly  solvable  Abel-type  integral  equations 
for  each  Fourier  coefficient  (a  la  A.  Cormack’s  work  on  x-ray 
CT).  While  the  2D  formulas  of  Ref.  [51]  involved  numeri¬ 
cal  instabilities,  this  was  fixed  in  the  3D  reconstructions  of 
Ref.  [52]. 

For  quite  some  time,  there  had  been  no  explicit  filtered 
backprojection-type  formula  obtained  for  any  closed  obser¬ 
vation  surface  S,  and  even  the  possibility  of  such  a  formula 
was  questioned.  Finally,  in  Ref.  [38],  such  formulas  were 
found  for  all  odd  dimensions  when  5  is  a  sphere,  and  then 
extended  to  even  dimensions  in  Ref  [53].  A  different  3D 
backprojection-type  formula  was  obtained  in  Ref  [54].  A 
backprojection  formula  for  arbitrary  dimension  was  found  in 
Ref  [55],  which  in  3D  coincided  with  the  one  in  Ref  [54].  It  is 
interesting  to  note  that  the  series  of  formulas  in  Refs.  [38,53] 


For  the  uniqueness  result  to  hold,  it  is  not  necessary  to  require  that  / 
vanishes  at  infinity,  sufficiently  fast  decay  at  infinity  also  suffices  145]. 
Ho\vever,  in  practical  situations, /does  vanish  at  large  distances. 


and  in  Refs.  [54,55]  are  not  equivalent  on  nonperfect  data, 
albeit  they  seem  to  work  numerically  equally  well  [30,31].  A 
different  derivation  of  formulas  in  Ref  [38]  has  recently  been 
proposed  [56].  No  closed  form  analytic  formulas  are  known 
for  surfaces  S  that  are  not  spheres,  cylinders,  or  planes. 

An  interesting  series  expansion  inversion  procedure  that 
theoretically  works  for  any  closed  observation  surface  S  was 
suggested  in  Ref.  [57].  The  expansion  of  the  image /(.r)  into 
the  eigenfunctions  of  the  Laplace  operator  A  inside  S  with 
zero  Dirichlet  conditions  on  S  is  obtained  in  terms  of  the  cor¬ 
responding  expansion  of  the  measured  data.  It  was  shown  in 
Ref  [57]  that  a  cubic  surface  used  as  S  rather  than  a  sphere 
works  much  better,  leading  to  very  fast,  accurate,  and  effi¬ 
cient  reconstructions.  This  series  expansion  procedure  was 
generalized  to  variable  sound  speeds  in  Ref.  [43],  albeit  it  is 
doubtful  that  this  can  lead  to  efficient  computations. 

Although  having  analytic  inversion  formulas  always  helps 
in  tomography,  it  is  also  known  that  reconstructions  can  often 
be  efficiently  done  without  having  explicit  formulas,  by  either 
algebraic  techniques,  or  their  combination  with  analytic  pre¬ 
conditioning.  This  is  true  in  TAT  as  well.  Namely,  it  is  not 
hard  to  write  a  good  approximate  inversion  operator  (tech¬ 
nically  called  a  parametrics)  that  gives  a  sufficiently  good 
approximate  reconstruction  and  preserves  the  locations  and 
strengths  of  all  singularities  of/ (e.g.,  sharp  boundaries).  Then 
one  can  bootstrap  the  quality  of  reconstruction  by  a  standard 
iterative  numerical  procedure  (e.g..  Refs.  [23,58,59]). 

Another  option  is  the  so-called  time  reversal.  Here  one 
solves  the  wave  equation  backward  in  time,  starting  with  a 
sufficiently  large  time,  when  the  signal  essentially  disappears. 
Then  when  time  t=0  is  reached,  the  image /(x)  is  recovered 
(e.g..  Refs.  [36-38,43,60,61]  and  references  therein). 

Various  other  discussions  of  inversion  procedures  and  rel¬ 
evant  references  can  be  found  in  Refs.  [29-31,36,37,62-71]. 

6.3  LIMITED  VIEW  PROBLEM 

We  now  switch  to  the  situation  of  our  main  interest  in  this 
paper — limited  view  TAT.  Suppose  that  S'  is  a  sphere  (or 
some  other  closed  surface)  surrounding  the  image  /  and  the 
only  possible  locations  of  transducers  belong  to  a  2D  piece 
T  of  S.  In  this  case,  we  only  have  the  data  g(y,t)  collected  at 
the  locations  y  in  T,  rather  than  from  the  whole  sphere.  In  this 
case,  we  face  a  limited  angle  (limited  view)  problem. 

6.3.1  Uniqueness  oe  Reconstruction 
EROM  Limited  View  Data 

Let  us  first  discuss  the  uniqueness  of  recovering  the  image 
from  the  data.  The  dimension  count  shows  that  one  should 
use  a  two-dimensional  piece  T  of  ^  in  order  to  expect  unique¬ 
ness.  Suppose  that  S  is  an  analytic  surface  (say,  a  sphere)  and 
T  its  two-dimensional  piece.  It  follows  from  the  results  of  Ref 
[35]  that  if  uniqueness  of  recovery  from  the  data  collected  on 
the  whole  S  is  known,  then  there  is  also  uniqueness  of  recov¬ 
ery  with  data  collected  on  T.  For  instance,  since  the  sphere  is 
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analytic  and  since  for  the  whole  sphere  uniqueness  of  recon¬ 
struction  is  known  (see  the  preceding  section),  we  conclude 
that  with  the  data  on  any  full-dimensional  piece  F  of  the 
sphere,  no  matter  how  tiny,  solution  of  the  TAT  problem  is 
unique.  The  same  holds  for  a  2D  piece  of  any  closed  analytic 
surface.^ 

Although  this  says  that  the  limited  view  data  is  theoreti¬ 
cally  sufficient,  anyone  trying  to  do  reconstructions  faces 
problems,  mostly  related  to  blurring  of  some  parts  of  the 
image.  This  effect  is  no  accident  and  can  be  predicted,  as 
explained  in  the  rest  of  the  chapter.  Let  us  look  first  at  how 
one  can  try  to  reconstruct  the  image  from  limited  view  data. 

6.3.2  Limited  View  Reconstructions 

The  simplest  thing  one  can  try  is  just  to  replace  the  missing 
data  with  zero  values* *  (the  procedure  often  called  zero-filling 
or  zero-padding)  and  then  reconstruct  the  image  as  if  the  data 
were  complete.  In  advance,  one  might  not  expect  any  reason¬ 
able  outcome  from  this,  however  the  following  observations 
appear,  independently  of  the  exact  procedure  involved: 

1.  Some  (we  will  call  them  “visible”^)  parts  of  the  tis¬ 
sue  interfaces  and  other  singularities  contained  in  the 
true  image  are  seen  clearly  and  at  the  right  locations. 
These  parts  do  not  depend  upon  the  exact  technique 
used,  but  do  depend  on  the  view  angle  available. 

2.  Other  parts  of  interfaces  blur  away. 

3.  “Smooth”  details  of  the  object,  i.e.,  exact  values  of 
f{x)  get  distorted,  in  some  areas  significantly. 

In  most  cases  in  tomography,  limited  data  problems  allow  no 
exact  reconstruction  formulas.  In  the  rare  cases  when  these 
do  exist,  this  does  not  change  the  effects  listed  above.  As  we 
explain  below,  there  is  a  good  reason  for  this,  and,  in  fact,  the 
“invisible”  and  “visible”  parts  of  singularities  can  be  easily 
predicted  by  a  simple  geometric  consideration.  In  order  to  do 
so,  we  need  to  introduce  some  technical  notions  from  the  so- 
called  microlocal  analysis  first  (see,  e.g..  Ref.  [72]  for  simple 
introduction  and  Refs.  [15,18,23,35,39,73,74]  for  applications 
in  integral  geometry  and  tomography). 

6.3.3  Wavefront  Sets  of  Functions 

Our  goal  is  to  apply  known  results  of  integral  geometry 
concerning  singularity  reconstruction  [11,12,15,39]  to  TAT 
[19,23].  The  exact  description  of  these  would  require  some 
notions  of  microlocal  analysis,  in  particular  the  notion  of  a 
wavefront  set  of  a  function  [72].  In  tomographic  problems, 
in  particular  in  TAT,  one  is  most  interested  in  only  one  type 


A  similar  result  holds  also  in  x-ray  CT,  where  limited  view  data  uniquely 
determine  the  image.  This  immediately  follows  from  the  projection-slice 
theorem  and  analyticity  of  the  Fourier  transform  of  a  compactly  sup¬ 
ported  function. 

*  It  is  advisable  sometimes  to  smooth  out  the  jump  between  zeros  and 
actual  data. 

*  Another  name  sometimes  used  is  “audible”  [39]. 


of  singularity:  the  jump  of/(x)  across  an  interface  (a  curve  in 
2D  or  a  surface  in  3D).  So,  we  will  describe  the  wavefront  set 
in  this  case  first.  Let/(x)  be  smooth  except  for  a  jump  across 
a  smooth  surface  L  (in  2D  case,  L  is  a  curve),  then  the  wave- 
front  set  WF(/)  of  f{x)  consists  of  pairs  (x,^),  where  point  x 
belongs  to  L  and  is  a  nonzero  vector  normal  to  L  at  x,  as 
shown  in  Figure  6.1. 

Before  introducing  the  wavefront  set  in  the  general  case,  let 
us  recall  that  smoothness  is  reflected  as  decay  in  the  Fourier 
domain.  Indeed,  if/(x)  is  smooth  and  has  compact  support  (or 
decays  sufficiently  fast  with  its  derivatives),  then  its  Fourier 
transform  /(^)  decays  faster  than  any  power  of  1^1  in  all 
directions  of  the  ^-space.  If  we  are  interested  in  detecting 
the  smoothness  of/ only  locally,  near  a  point  Xq,  we  cut  other 
parts  off  by  multiplying /by  a  smooth  function  (|)  that  is  non¬ 
zero  only  near  Xq.  Then,  again,  the  Fourier  transform  (|)/(^) 
of  the  product  decays  faster  than  any  power  of  1^1  in  all  direc¬ 
tions  of  the  ^-space.  Well,  what  does  it  mean  now  that/is  not 
smooth  near  Xq?  This  means  that  (j)/(^)  will  not  be  decaying, 
unless  (|)  vanishes  at  Xq,  which  we  will  prohibit.  However,  it 
might  still  decay  in  some  directions.  This  now  leads  to  the 
general  definition  (e.g..  Refs.  [39,72])  of  the  wavefront  set, 
which  picks  up  for  each  point  Xg  the  bad  directions  ^  only.  It 
thus  consists  of  pairs  (x,^),  where  ^  must  be  a  nonzero  vec¬ 
tor  (in  order  to  determine  a  direction).  It  is  easier  to  describe 
which  pairs  (xg,^)  (where  ^^0)  do  not  belong  to  the  wave- 
front  set  WF(/)  of  a  function/  Namely,  this  happens  if  there 
is  a  neighborhood  u  of  Xg  and  a  conic  neighborhood 


V  such  that 


A  iiL 

1^1  l^ol 


of  the  direction  such  that  for  any  smooth  function  (|)  sup¬ 
ported  inside  u,  the  Fourier  transform  (|)/(^)  of  the  function 
([/decays  to  zero  faster  than  any  power  of  1^1  when 
in  V^.  The  role  of  the  function  (j)  is  to  eliminate  the  possi¬ 
ble  bad  behavior  of  /  away  from  Xg  and  keep  only  the  local 


FIGURE  6.1  If/(x)  is  smooth,  except  a  jump  singularity  across  L, 
then  its  wavefront  set  WF(/)  consists  of  pairs  (x,^),  where  x  belongs 
to  the  jump  interface  L,  and  ^  is  a  nonzero  vector  normal  to  L,  at 
X.  (Reproduced  from  Xu,  Y.,  L.V.  Wang,  G.  Ambartsoumian,  and 
R  Kuchment,  Med.  Phys.  31(4):724-33,  2005.  With  permission.) 
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information  about  properties  of  /  at  Xq.  Thus,  the  wavefront 
set  of/ carries  the  information  about  singularities  of/  which 
is  localized  hoth  in  spatial  variable  x  and  frequency  variable 
^  (the  word  used  is  “microlocalization”,  and  thus  “microlocal 
analysis”).  Projecting  the  wavefront  set  onto  the  space  com¬ 
ponent  X  (by  forgetting  the  frequency  variable  ^),  one  obtains 
the  so-called  singular  support  of/  i.e.,  the  set  of  all  singulari¬ 
ties  of/ — points  near  which /is  not  smooth.  For  instance,  it 
is  easy  to  check  that  the  wavefront  set  of  the  delta  function 
consists  of  (0,  for  any  ^^0. 

6.3.4  "Visible"  Singularities  and  Their 
Reconstruction 

The  general  idea  of  the  microlocal  approach  to  limited  view 
problems  is  the  following  [11]  (see  a  somewhat  different,  more 
limited  approach  in  Ref.  [12]).  One  tries  to  determine  which 
wavefront  set  elements  (xfy  of  the  object /lead  to  singulari¬ 
ties,  i.e.,  wavefront  set  elements  in  the  measured  data  g.  This 
can  be  done  by  the  technique  of  the  so-called  Fourier  integral 
operators  (FIOs)  [73,74],  which  is  beyond  the  scope  of  the 
current  consideration.  If  such  a  pair  (xfy  does  lead  to  a  singu¬ 
larity  in  g,  it  is  called  “visible”  or  “audible”.  Such  wavefront 
set  points  of  the  object  can  be  stably  reconstructed  from  the 
data  g,  and  thus  the  singularity  will  show  up  in  the  recon¬ 
struction.  If  a  particular  wavefront  set  point  does  not  influ¬ 
ence  the  singularities  of  the  data  g  (i.e.,  it  is  smoothed  out),  it 
becomes  “invisible”  and  thus  will  be  blurred  away.  This  blur¬ 
ring  effect  is  intrinsic  to  the  problem  and  cannot  be  overcome 
by  changing  analytic  or  numerical  techniques,  unless  some 
extra  information  is  incorporated  into  the  problem. 

This  concept  is  not  hard  to  understand.  Imagine,  for 
instance,  that  the  operator  that  maps  the  unknown  /  into 
the  data  set  g,  turns  all  images  even  the  ones  with  singu¬ 
larities,  into  smooth  functions  g.  In  the  Fourier  domain,  this 
can  be  interpreted  as  an  operation  that  turns  slowly  decaying 
Fourier  transforms  into  the  ones  that  decay  faster  than  any 
power.  For  instance,  one  can  imagine  that  this  is  done  by 
applying  a  filter  that  decays  extremely  fast,  and  essentially 
acts  as  a  low-pass  filter.  Trying  to  invert  the  procedure  and 
reconstruct  the  object/from  the  data  g,  one  runs  into  trouble, 
since  the  required  filters  will  grow  faster  than  any  power 
(often  exponentially).  This  clearly  indicates  impossibility 
of  stable  reconstruction  of  high-frequency  components,  and 
thus  of  any  sharp  details.^  A  “microlocal”  extension  of  this 
consideration  shows  that  if  some  part  of  the  wavefront  set 
(i.e.,  a  singularity)  of  the  image /does  not  contribute  to  the 
wavefront  set  of  the  data,  then  this  singularity  cannot  be  sta¬ 
bly  recovered  from  the  data. 

It  is  thus  important  to  be  able  to  understand  in  advance 
which  wavefront  set  elements  (x,^)  of  the  (unknown) /would 
lead  to  some  wavefront  set  elements  in  g.  This  would  deter¬ 
mine  which  singularities  are  “visihle”  from  the  data  g. 


^  One  faces  such  instabilities,  for  instance,  trying  to  solve  the  heat  equation 
back  in  time,  or  dealing  with  electrical  impedance  imaging  or  optical 
tomography. 


Fortunately,  there  is  a  very  simple  answer  to  this  question 
in  the  cases  of  x-ray  tomography,  SPECT,  TAT,  and  some 
other  imaging  methods  [11,15,18,19,39]  (while  its  justification 
is  very  nontrivial).  We  will  describe  it  following  Ref.  [19]  (see 
also  Ref.  [23]),  as  applied  to  TAT: 

The  “visibility”  condition 

An  element  (x,  ^)  of  the  wavefront  set  of  f  cannot  be  recov¬ 
ered  looking  at  the  singularities  of  the  spherical  integrals 
data,  unless  there  is  a  detector  location  and  a  sphere  cen¬ 
tered  at  this  location  that  passes  through  the  point  x  and  is 
normal  to  ^  at  this  point.  In  other  words,  in  TAT  one  can  see 
without  blurring  only  those  parts  of  the  interfaces  that  can 
be  touched  tangentially  by  spheres  centered  at  detector  posi¬ 
tions.  This  means  that  in  order  to  recover  stably  (i.e.,  without 
blurring)  the  interface  L,  one  needs  to  have  for  each  point  of 
L  a  detector  located  along  the  normal  line  to  L  at  this  point. 

If  for  some  part  of  L,  the  normals  do  not  pass  through  the 
detectors,  this  part  will  be  “invisible”  and  will  be  mandato- 
rily  blurred  away. 

Remark  1  This  principle  does  not  depend  upon  a  particular 
reconstruction  method.  So,  a  bad  method  can  increase  blur¬ 
ring,  but  even  the  best  methods  cannot  recover  sharply  the 
“invisible”  interfaces.  This  is  why  in  this  text  we  do  not  go 
into  any  details  of  reconstructions.  Certainly,  incorporation 
of  additional  a  priori  information  about  the  image  (e.g.,  that 
it  consists  of  a  few  simple  “blocks”)  could  potentially  lead 
to  improvements. 

The  visibility  condition  described  above  is  not  hard  to 
understand.  Indeed,  assume  for  simplicity  that  the  integration 
that  produces  the  measured  data  is  done  along  planes  rather 
than  spheres,  and  that  the  interface  is  also  planar.  Then,  if 
one  integrates  along  a  stack  of  parallel  planes  x  03=s  that  is 
not  parallel  to  the  interface  L  (i.e.,  the  normal  vector  co  to  the 
plane  of  integration  is  not  normal  to  L),  the  singularity  of  / 
along  the  interface  is  smoothed  out,  and  the  resulting  value 
depends  on  s  smoothly.  Only  if  one  has  at  one’s  disposal  a  set 
of  planes  parallel  to  the  interface  (i.e.,  co  is  normal  to  L),  then 
the  integration  of/ will  result  in  a  singularity  with  respect  to 
s.  This  indicates  that  invisible  parts  of  the  wavefront  set  do 
not  show  up  as  singularities  in  the  measured  data,  and  thus 
cannot  be  stably  reconstructed.  This  hand-waving  explana¬ 
tion  can  be  made  precise,  with  substantial  technical  effort. 

Notice  that  if  at  some  location  x,  any  line  passing  through 
X  crosses  a  detector  position,  then  for  any  image,  we  expect 
its  details  near  x  to  be  reconstructed  sharply.  This  leads  us  to 
the  following  notion: 

Definition  2  Suppose  that  detectors  can  be  located  along  a 
2D  (ID  in  the  planar  case)  piece  F  of  the  observation  surface 
S  only.  The  detectable  region  consists  of  such  points  x  inside 
S  that  every  line  passing  through  x  crosses  F 

According  to  the  visibility  condition,  any  object  supported 
inside  the  detectable  region  will  be  sharply  reconstructed, 
while  the  objects  reaching  outside  the  detectable  region  will 
have  some  parts  of  their  interfaces  blurred  away. 
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FIGURE  6.2  (a)  Illustration  of  the  “detectable  regions”  (shaded  areas)  of  circular  Radon  transformation,  when  the  detector  moves  along  a 

single  arc  (solid)  of  a  circle,  (b)  Two  arcs,  (c)  Three  arcs.  (Reproduced  from  Xu,  Y.,  L.V.  Wang,  G.  Ambartsoumian,  and  R  Kuchment,  Med. 
Phys.  31(4):724-33,  2005.  With  permission.) 


FIGURE  6.3  (a)  “Visible”  (solid  line)  and  “invisible”  (dashed)  boundaries  of  a  square  object,  and  the  “detectable  regions”  (shaded  areas) 

when  the  detector  moves  along  an  arc  (solid),  (b)  Same  as  (a)  for  a  disk  phantom,  (c)  Same  as  (a)  except  that  the  detector  moves  along  the  line 
segment  AB  and  the  objects  are  a  square  and  a  disk.  The  “visible”  boundaries  are  expected  to  be  recoverable  stably,  while  the  “invisible” 
boundaries  should  be  blurred  away.  (Reproduced  from  Xu,  Y.,  L.V.  Wang,  G.  Ambartsoumian,  and  P.  Kuchment,  Med.  Phys.  31(4):724-33, 
2005.  With  permission.) 


Figure  6.2  illustrates  the  detectable  regions.  Figure  6.3 
shows  the  expected  behavior  of  the  reconstructed  interfaces 
of  simple  square  and  disk  phantoms  with  limited  view  data. 

Results  of  reconstructions  from  synthetic  and  experi¬ 
mental  data,  which  confirm  this  conclusion,  are  shown  in 
Section  6.4. 


gested:  better  approximate  inverses,  corrective  coefficients, 
numerical  minimization,  using  range  conditions  for  recover¬ 
ing  the  missing  data,  etc.  (e.g..  Refs.  [22,26,27]).  A  very  recent 
work  [77]  shows  great  promise  for  the  final  resolution,  albeit 
at  the  moment  of  writing  this  chapter,  not  all  necessary  details 
have  been  filled  in. 


6.3.5  Half-Sphere  Problem 

A  particular  case  of  interest,  which  has  attracted  the  attention 
of  several  researchers,  is  when  the  detectors  can  be  placed 
along  half  of  the  observation  sphere  S,  and  thus  the  detec¬ 
tion  region  is  the  corresponding  half-ball.  It  is  assumed  that 
the  object  is  located  inside  the  detection  region.  Then,  the 
principle  discussed  in  the  previous  section  predicts  that  all 
singularities  of  the  object  are  “visible”  and  thus  reconstruc¬ 
tions  should  be  sharp.  If  one  attempts  to  zero-fill  the  data 
from  the  other  half  sphere  and  use  any  of  the  standard  recon¬ 
struction  methods,  the  interfaces  truly  remain  sharp,  but  the 
intensities  deteriorate  towards  the  missing  equator.  Having  an 
exact  reconstruction  formula  for  half-sphere  data  would  fix 
this  problem,  but  such  a  formula  has  not  been  found  so  far.” 
So,  different  partial  remedies  for  this  ailment  have  been  sug- 


In  SPECT,  a  similar  problem  waited  quite  a  long  time  before  its  satisfac¬ 
tory  resolution  175,76]. 


6.3.6  Local  Tomography  and 
Singularity  Sharpening 

We  would  like  to  mention  briefly  the  principle  of  the  so- 
called  local  tomography  [18,78-80].  In  this  method,  before 
backprojection,  an  additional  growing  filter  in  the  frequency 
domain  is  applied  in  order  to  sharpen  singularities.  The  result¬ 
ing  reconstruction  has  incorrect  numerical  values  of/(x),  but 
significantly  emphasized  interfaces  and  other  singularities, 
which,  for  instance,  can  be  useful  when  small  blood  vessels 
or  region  of  interest  tomography  are  of  interest.  Local  tomog¬ 
raphy  applies  to  TAT  as  well  [19,23].  In  the  case  of  limited 
view  data,  it  also  recovers  the  “visible”  parts  of  the  interfaces 
only.  Some  of  the  reconstructions  shown  in  the  next  section 
include  their  local  tomography  counterparts. 

It  is  interesting  to  notice  (see  more  about  this  in  the  fol¬ 
lowing  section),  that  transducers’  responses  often  act  as  that 
extra  filter  needed  for  local  tomography,  and  thus  boundaries 
are  emphasized  without  any  extra  effort. 
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6.4  NUMERICAL  AND  EXPERIMENTAL 
VERIFICATION 

In  this  section,  we  will  illustrate  our  theoretical  analysis  and 
conclusions  with  reconstructions  [23]  from  both  synthetic  and 
experimental  (i.e.,  collected  from  a  physical  phantom)  TAT 
data.  We  do  not  mention  details  of  specific  reconstruction 
methods  used  because,  as  explained  earlier,  the  “visibility- 
invisibility”  effect  does  not  depend  on  the  method  used. 

6.4.1  Reconstruction  from  Synthetic  Limited 
View  Thermoacoustic  Tomography  Data 

A  numerical  phantom  that  contains  four  sharp  and  one  soft 
inclusion  is  shown  in  Figure  6.4.  Among  the  sharp  inclusions, 
we  have  one  large  and  two  small  squares  and  one  disk.  The 
object  value  is  set  to  be  0.5  within  the  largest  square,  unity 
within  other  sharp  inclusions,  and  zero  elsewhere.  Inside  the 
“soft”  circular  inclusion,  this  value  drops  linearly  with  the 
radius  from  unity  at  the  center  to  zero  at  the  interface,  in  order 
to  simulate  a  gradual  interface.  The  imaged  field  of  154  mm  by 
154  mm  is  mapped  with  a  128  x  128  mesh.  The  detection  circle 
has  a  radius  of  133  mm  and  is  centered  at  the  center  of  the  pic¬ 
ture.  We  scan  200  steps  in  all  the  simulations.  The  gray  scale 
and  the  scale  bar  of  the  images  are  shown  below  the  images  in 
Figure  6.5.  The  top  row  of  reconstructions  employs  the  local 
tomography  formula  that  emphasizes  the  boundaries.  The  next 
one  uses  the  approximate  filtered-backprojection  (FBP)  for¬ 
mula  in  Ref  [23],  and  the  lowest  one  shows  the  improvements 


FIGURE  6.4  Diagram  of  inclusions  in  TAT  (used  in  Figure  6.5). 
The  value  of  the  image /(x)  is  set  to  he  0.5  in  the  largest  square,  and 
unity  within  other  sharp  inclusions,  and  zero  elsewhere.  Inside  the 
“soft”  circular  inclusion,  this  value  drops  linearly  with  the  radius 
from  unity  at  the  center  to  zero  at  the  interface.  (Reproduced  from 
Xu,  Y.,  L.V.  Wang,  G.  Amhartsoumian,  and  R  Kuchment,  Med. 
Phys.  31(4):724-33,  2005.  With  permission.) 


achieved  by  running  the  algebraic  reconstruction  method 
(TCG),  starting  with  the  FBP  as  an  initial  guess  [23]. 

The  left  column  uses  only  the  data  collected  from  the 
Till  detection  arc  in  the  first  quadrant.  None  of  the  phantom 
inclusions  fit  into  the  “detectable  region”.  One  can  see  that 
all  parts  of  the  inclusion  boundaries  the  normals  to  which  do 
not  intersect  the  detection  arc  are  blurred  (even  in  the  local 
tomography  reconstruction).  Other  parts  of  the  boundaries 
are  sharp.  This  is  in  perfect  agreement  with  our  theoretical 
prediction.  The  soft  inclusion  is  not  significantly  affected  by 
the  artifacts. 

The  middle  column  employs  the  data  collected  from  the 
detection  arc  of  approximately  217  degrees  (the  angle  0  in 
Figure  6.4),  whose  chord  coincides  with  the  bottom  side 
of  the  large  square  inclusion.  In  this  case,  all  inclusions 
are  in  the  “detectable  region”,  and  hence  all  the  boundar¬ 
ies  are  reconstructed  sharply.  The  third  column  represents 
the  full  data  reconstruction.  Notice  that  the  quality  of  the 
final  reconstructions  in  the  last  two  columns  is  the  same. 
Figure  6.6a  and  b  show  the  reconstructed  image /(x)  along 
the  dashed-dotted  line  in  Figure  6.4,  using  the  FBP  (Figures 
6.5d  through  f)  and  TCG  reconstructions  (Figure  6.5g 
through  i),  respectively.  The  exact  value  is  also  shown  for 
comparison.  It  can  be  found  in  Figure  6.6a  that  the  results 
of  FBP  are  in  good  agreement  with  the  real  value  for  the 
case  of  217-degree  and  360-degree  detection,  where  all 
objects  are  in  the  “detectable  region”.  Iteration  improves  the 
results  further,  as  shown  in  Figure  6.6b.  Even  for  the  case 
of  a  90-degree  detection  curve,  the  profile  of  the  objects  is 
reconstructed.  Comparing  Figure  6.6a  and  b,  one  finds  that 
the  significant  overshoot  and  undershoot  in  FBP  can  be  con¬ 
siderably  reduced  by  TCG  iterations  (we  remind  the  reader 
that  FBP  is  only  an  approximation  rather  than  the  implemen¬ 
tation  of  an  exact  formula). 

Figure  6.7  shows  the  relative  error  of  each  reconstruc¬ 
tion  as  a  function  of  the  scanned  angular  range  with  respect 
to  the  center  of  the  scan.  We  study  the  mean  reconstruction 
values  in  the  hard  sphere,  the  central  square,  and  the  back¬ 
ground.  The  errors  of  reconstruction  are  normalized  to  the 
corresponding  real  values  in  the  cases  of  the  hard  sphere 
and  the  central  square,  and  to  the  real  value  of  the  hard 
sphere  in  the  case  of  the  background  (because  its  real  value 
is  zero).  When  the  scanned  angular  range  is  less  than  n, 
the  errors  decrease  sharply  with  increasing  scanned  angu¬ 
lar  range.  By  contrast,  when  the  scanned  angular  range  is 
larger  than  Ji,  the  errors  change  much  more  slowly  as  the 
scanned  angular  range  increases.  The  results  agree  with 
our  theoretical  conclusions.  However,  there  are  some  fluc¬ 
tuations  added  to  the  trends  of  the  curves.  By  comparing 
the  three  curves  in  Figure  6.7,  we  find  that  these  fluctua¬ 
tions  depend  strongly  on  the  location  of  the  object  with 
respect  to  the  detection  curve.  More  extensive  study  is 
needed  to  understand  these  fluctuations.  There  are  some 
residual  errors  even  in  the  full-view  detection  in  Figure  6.7. 
This  is  because  we  used  an  approximate  backprojection 
algorithm,  rather  than  an  exact  inversion  (which  was  not 
available  at  that  time). 
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FIGURE  6.5  Images  reconstructed  from  simulated  TAT  data  corresponding  to  the  phantom  in  Figure  6.4.  The  three  columns  correspond 
from  left  to  right  to  detection  angles  of  90  degrees  (from  0°  to  90°),  217  degrees  (from  -19°  to  198°  as  shown  hy  the  angle  0  in  Figure  6.4), 
and  360  degrees,  respectively.  The  three  rows  correspond  from  top  to  bottom  to  the  local  tomographic  reconstruction,  FBP,  and  FBP  with 
the  consecutive  TCG,  respectively.  The  values  (minimum,  maximum)  of  the  gray  scale  for  (a-i)  are  (-0.8081,  1.0000),  (-0.8302,  1.0000), 
(-0.7515,  1.0000),  (-2.0745,  1.7899),  (-0.6385,  1.0723),  (-0.1030,  1.0349),  (-0.9284,  1.2859),  (-0.0326,  1.0030),  and  (-0.0149,  1.0021), 
respectively.  The  maxima  of  the  local  reconstructions  are  normalized  to  unity.  (Reproduced  from  Xu,  Y.,  L.V.  Wang,  G.  Ambartsoumian, 
and  P.  Kuchment,  Med.  Phys.  31(4):724-33,  2005.  With  permission.) 


6.4.2  Reconstruction  from  Experimental  Limited 
View  Thermoacoustic  Tomography  Data 

The  experimental  setup  is  described  in  [23,67]  and  will 
not  be  repeated  here.  The  sample  and  the  polar  coordinate 
system  describing  the  scanning  orbit  are  shown  in  Figure 
6.8a.  The  sample  consists  of  a  muscle  cylinder  of  4  mm  in 
diameter  and  5  mm  in  length,  embedded  in  a  chunk  of  pork 


fat  of  1.2  cm  in  radius,  There  is  a  10-mm  fat  layer  below 
the  muscle  and  another  7-mm  one  above  it.  An  EM  pulse 
is  delivered  to  the  sample  from  below  (i.e.,  from  behind 
the  picture  plane).  With  a  scanning  radius  of  r^=1.l  cm, 
thermoacoustic  data  are  collected  around  the  sample  over 
a  2n  angular  span  with  161  steps.  The  EM  pulse  profile 
and  the  impulse  response  function  of  the  ultrasonic  trans¬ 
ducer  impose  a  filter  on  the  thermoacoustic  signals.  We 
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FIGURE  6.6  (a)  The  graphs  of  FBP  reconstructions  shown  in  Figure  6.5d  through  f,  and  the  corresponding  exact  value  along  the  dashed- 

dotted  line  in  Figure  6.4.  (b)  The  graphs  corresponding  to  TCG  reconstructions,  Figure  6.5d  through  f,  along  the  same  line  as  in  (a). 
(Reproduced  from  Xu,  Y.,  L.V.  Wang,  G.  Ambartsoumian,  and  R  Kuchment,  Med.  Phys.  31(4):724-33,  2005.  With  permission.) 


FIGURE  6.7  Dependence  of  the  relative  errors  of  the  mean  values 
in  the  hard  sphere  (circle  markers),  the  central  square  (square  mark¬ 
ers),  and  the  background  (asterisks)  on  the  scanned  angular  range. 
(Reproduced  from  Xu,  Y.,  L.V.  Wang,  G.  Ambartsoumian,  and  P. 
Kuchment,  Med.  Phys.  31(4):724-33,  2005.  With  permission.) 

attempted  to  correct  this  effect  using  deconvolution,  but 
found  that  the  resulting  images  were  distorted,  due  to  the 
lack  of  precise  knowledge  of  the  filter.  Therefore,  we  do 
not  use  deconvolution  in  the  reconstruction.  This  leads, 
as  explained  earlier,  to  somewhat  emphasized  interfaces. 
Figure  6.8b  through  d  show  the  reconstructed  images  using 
FBP  with  three  sets  of  data.  In  Figure  6.8b,  we  choose  the 
data  collected  along  a  circular  detection  arc  of  92  degrees, 
located  at  the  top  of  the  picture  and  almost  symmetric  with 
respect  to  its  vertical  axes.  One  sees  that  the  left  and  right 


boundaries  of  the  muscle  cylinder  and  of  the  pork  chunk 
are  blurred  away,  since  their  normal  lines  do  not  touch  the 
detection  arc,  while  the  rest  of  the  boundary  is  sharp.  The 
next  figure  shows  the  reconstructed  image  obtained  with 
the  data  collected  from  a  202-degree  arc,  where  the  whole 
phantom  fits  into  the  detectable  region.  All  boundaries  are 
sharp  now.  Finally,  the  last  figure  shows  the  image  recon¬ 
structed  with  the  full-view  data. 

Notice  that  although  no  local  reconstruction  algorithms 
are  applied,  the  boundaries  are  somewhat  emphasized.  The 
reason  for  this  is  the  presence  in  the  data  of  the  impulse 
response  function  of  the  ultrasonic  transducer,  which  has  an 
effect  similar  to  the  application  of  an  additional  derivative 
with  respect  to  the  radius  of  the  circle  of  integration.  The 
presence  of  such  a  derivative  emphasizes  high  frequencies 
and  makes  the  reconstruction  similar  to  a  version  of  a  local 
tomography  algorithm. 

6.4.3  Discussion  on  Experimental  Results 

As  mentioned  earlier,  although  circular  scanning  is  used  in 
both  our  numerical  and  experimental  studies,  our  conclusions 
can  be  applied  to  other  configurations  as  well.  In  TAT  with  a 
planar  configuration  [62,64-66],  detections  are  implemented 
on  a  part  of  a  line  or  a  plane  where  the  scanning  view  is  quite 
limited;  consequently,  artifacts  and  interface  blurring  appear 
in  the  reconstructed  images.  In  fact,  in  planar  and  linear 
scanning  geometries,  one  can  never  have  an  object  immersed 
entirely  into  the  “detectable  region”  because  the  normal  lines 
to  any  interfaces  that  are  orthogonal  to  the  detection  plane 
(line)  never  pass  through  a  detector.  As  a  consequence,  those 
parts  of  the  interfaces  will  be  blurred  in  any  kind  of  recon¬ 
struction.  For  a  sufficiently  large  view,  these  parts  will  be 
small,  but  theoretically  will  never  disappear.  For  example,  2D 
planar  detection  is  utilized  to  image  artificial  blood  vessels 
[64];  the  scanning  view  is  about  2.18  steradians.  Therefore,  it 
is  not  surprising  that  only  the  interfaces  more-or-less  parallel 
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FIGURE  6.8  (a)  Photograph  of  the  experimental  sample,  (b-d)  TAT  images  reconstructed  using  detection  arcs  of  92  degrees  (from  50°  to 

142°  in  (a)),  202  degrees  (from  -18°  to  184°),  and  360  degrees,  respectively.  The  blurred  parts  of  the  boundaries  in  (b)  due  to  the  limited  view 
agree  with  the  theoretical  predictions.  In  (c)  all  the  boundaries  are  resolved,  since  the  object  fits  into  the  “detectable  region”.  (Reproduced 
from  Xu,  Y.,  L.V.  Wang,  G.  Ambartsoumian,  and  P.  Kuchment,  Med.  Phys.  31(4):724-33,  2005.  With  permission.) 


to  the  plane  of  detection  are  well  imaged.  Linear  scanning 
detection  was  used  in  Ref.  [66]  to  image  a  2D  phantom. 
Because  the  view  of  the  linear  scanning  in  Ref.  [66]  is  much 
larger  than  that  of  planar  scanning  in  Ref.  [64],  the  inter¬ 
faces  are  recovered  much  more  completely.  However,  due  to  a 
limited  view,  artifacts  and  interface  blurring  similar  to  those 
demonstrated  in  our  numerical  and  experimental  studies  still 
appear  in  the  images  [66]. 

By  comparing  Figures  6.5  and  6.8,  we  observe  that  the 
images  reconstructed  from  incomplete  data  when  an  object 
is  in  the  detectable  region,  have  comparable  quality  with 
those  from  the  full-view  data.  Scanning  a  smaller  range  has 
the  advantages  of  reducing  the  scanning  time  or  the  size  of 
the  acoustic  transducer  array.  It  should  be  pointed  out  that 
this  advantage  usually  exists  when  both  the  sample  and 
medium  are  relatively  acoustically  homogeneous.  When 


strong  wavefront  distortion  caused  by  acoustic  heterogene¬ 
ities  occurs,  it  might  be  beneficial  to  collect  signals  from  all 
directions. 

6.5  ADDITIONAL  REMARKS  AND 
CONCLUSIONS 

As  mentioned  before,  one  may  incorporate  some  additional 
information  about  the  image,  or  change  the  physical  set-up 
of  the  problem  to  stabilize  the  inverse  problem  and  make  all 
or  some  formerly  invisible  interfaces  visible.  Recently,  it  was 
shown  [70]  that  taking  into  account  some  a  priori  knowledge 
about  the  interfaces  leads  to  reconstruction  of  previously 
invisible  parts.  In  another  direction,  acoustic  reflectors  were 
proposed  as  a  means  of  reflecting  the  acoustic  waves,  which 
would  otherwise  not  be  measured,  back  onto  the  sensor.  It  was 
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shown  in  that  case  an  existing  FFT-based  image  reconstruc¬ 
tion  algorithm  can  be  used  to  reconstruct  the  image  without 
the  limited  view-induced  blurring  [71,81]. 

The  main  points  of  this  survey  can  be  summarized  as 
follows: 

•  A  geometric  principle  is  described  that  allows  a  sim¬ 
ple  determination  of  which  sharp  parts  of  the  object 
are  expected  to  be  blurred  when  reconstructed  from 
limited  view  thermoacoustic  data. 

•  This  blurring  is  independent  of  the  particular  recon¬ 
struction  method  and  cannot  be  overcome,  unless 
some  extra  information  about  the  object  is  known. 

•  Numerical  results  using  synthetic  and  experimental 
data  are  shown  that  support  the  conclusions. 
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